load '../data/DataforMatlab/UK_Budget_S'
load '../data/DataforMatlab/pvec_tmp.mat'

pvec_tmp=pvec_tmp(end-16:end,:);
B_vec = B_vec./sum(B_vec);
I = size(B_vec,1);N = size(B_vec,2);T = size(B_vec,3);
pvec = reshape(pvec_tmp',[I,1,T]);
[U_vec] = CalMoneyMetric(I_vec, B_vec, pvec,0);
 

load ../data/DataforMatlab/IP.csv
load ../data/DataforMatlab/torn_CPI.csv

IP = IP(:,end-16:end);
if exist('PC', 'var') == 1
    PC = PC(:,end-16:end);
end
%% Percentile for plotting
% Reshape data
Wfull = reshape(I_vec(:,:,:), [N,T]);
U = reshape(U_vec,[N,T]);
% Initialize variables
PC = zeros(9,T);
Wq = zeros(9,T);
Uq = zeros(9,T);

for p = 1:9
    PC(p,:) = sum(Wfull <= IP(p,:));
    PC(PC == 0) = 1;
    Wq(p,:) = Wfull(sub2ind(size(Wfull), PC(p,:), 1:T));
    Uq(p,:) = U(sub2ind(size(U), PC(p,:), 1:T));
end

%%Insert NAN into income vectors for which the money metric cannot be calculated.
[~, indmax] = max(U);
[~, indmin] = min(U);
W = NaN(size(Wfull));
Wmax = Wfull(sub2ind(size(Wfull), indmax, 1:T));
Wmin = Wfull(sub2ind(size(Wfull), indmin, 1:T));
for tau = 1:T
    W(indmin(tau):indmax(tau), tau) = Wfull(indmin(tau):indmax(tau), tau);
end

logUmax = max(log(U));
logUmin = min(log(U));

RC = [Wmin;Wq;Wmax]./exp(torn_CPI') * 52;
OurxaxisUK = [Wmin(:,end);Wq(:,end);Wmax(:,end)]*52;
OuryaxisUK = [exp(logUmin(:,end)+log(52)) ;(Uq(:,end)*52);exp(logUmax(:,end)+log(52)) ];



WelfareRelevant = log( [min(W)./min(U);Wq./Uq;max(W)./max(U)]);
diffp = WelfareRelevant - torn_CPI';

agg_OuryaxisUK= OuryaxisUK;
agg_OurxaxisUK = OurxaxisUK;
agg_U_vec = U_vec;
agg_I_vec = I_vec;
agg_diffp = diffp;
save("../data/DataforMatlab/agg_mmcpi.mat",'agg_OuryaxisUK','agg_OurxaxisUK','agg_U_vec','agg_I_vec','agg_diffp');

    function F_X = interpNaN(X,Y)
        indexNaN = isnan(X);
        X(indexNaN) = [];% Eliminate NaN
        Y(indexNaN) = [];% Eliminate NaN
        [X, indexUN] = unique(X); 
        F_X = griddedInterpolant(X,Y(indexUN));
    end
